Disordered bosons in one dimension: from weak to strong randomness criticality 
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We investigate the superfluid-insulator quantum phase transition of one-dimensional bosons with 
off-diagonal disorder by means of large-scale Monte-Carlo simulations. For weak disorder, we find 
the transition to be in the same universality class as the superfluid-Mott insulator transition of 
the clean system. The nature of the transition changes for stronger disorder. Beyond a critical 
disorder strength, we find nonuniversal, disorder-dependent critical behavior. We compare our 
results to recent perturbative and strong-disorder renormalization group predictions. We also discuss 
experimental implications as well as extensions of our results to other systems. 



Bosonic many-particle systems can undergo quan- 
tum phase transitions between superfluid and localized 
ground states due to interactions and lattice effects. 
These superfluid-insulator transitions occur in a wide 
variety of experimental systems ranging from helium in 
porous media, Josephson junction arrays, and granular 
superconductors to ultracold atomic gases [iHll . In many 
of these applications, the bosons are subject to quenched 
disorder or randomness. Understanding the effects of dis- 
order on the superfluid-insulator transition and on the 
resulting insulating phases is thus a prime question. 

The case of one space dimension is especially inter- 
esting because the superfluid phase is rather subtle and 
displays quasi- long-range order instead of true long-range 
order. Moreover, the Anderson localization scenario for 
non-interacting bosons suggests that disorder becomes 
more important with decreasing dimensionality. 

Giarmarchi and Schulz Q studied the influence of weak 
disorder on the interacting superfluid by means of a per- 
turbative renormalization group analysis. They found 
the superfluid-insulator transition to be of Kosterlitz- 
Thouless (KT) type 0, with universal critical expo- 
nents and a universal value of the Luttinger parameter 
g = tt^/PsK at criticality (p s denotes the superfluid stiff- 
ness and k the compressibility). This analysis was re- 
cently extended to second order in the disorder strength, 
with unchanged conclusion [111 ]. 

A different scenario emerges, however, from the real- 
space strong-disorder renormalization group approach. 
In a series of papers [12J, Altman et al. studied one- 
dimensional interacting lattice bosons in various types 
of disorder. In all cases, they found that the superfluid- 
insulator transition is characterized by KT-like scaling 
of lengths and times, but it occurs at a nonuniversal, 
disorder-dependent value of the Luttinger parameter. 
The transition is thus in a different universality class than 
the weak-disorder transition 0. However, Monte-Carlo 



simulations 13] did not find any evidence in favor of the 
strong-disorder critical point. 

In view of these seemingly incompatible results, it is 
important to determine whether or not both types of 
superfluid-insulator critical points indeed exist in sys- 
tems of interacting disordered bosons in one dimension. 
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FIG. 1. (Color online) Critical Luttinger parameter g and ex- 
ponent rj [plotted as l/(2r])] of the superfluid-insulator tran- 
sition as functions of the disorder strength 1 — r. The critical 
behavior appears universal for weak disorder while it becomes 
disorder-dependent for strong disorder. The lines are guides 
to the eye only. 



Moreover, it is important to study whether they can be 
reached for realistic disorder strengths. 

In this Letter, we employ large-scale Monte-Carlo sim- 
ulations to address these questions. We focus on the 
case of off-diagonal disorder at large commensurate fill- 
ing; other types of disorder will be discussed in the con- 
clusions. Our results can be summarized as follows (see 
Fig.[T]). For weak disorder, we find a KT critical point in 
the universality class of the clean (1+1) dimensional XY 
model, with universal exponents and a universal value of 
the Luttinger parameter at the transition. This agrees 
with the predictions of the perturbative renormalization 
group. If the disorder strength is increased beyond a 
threshold value, the nature of the transition changes. 
While the scaling of length and time scales remains KT- 
like, the critical exponents and the Luttinger parame- 
ter become non-universal, in agreement with the strong- 
disorder scenario 



12j. In the remainder of this Letter, 



we explain how these results were obtained, we discuss 
their generality as well as implications for experiment. 

The starting point is the disordered one-dimensional 
quantum rotor Hamiltonian 
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which represents, e.g., a chain of superfiuid grains with 
Josephson couplings Jj, charging energies Uj and offset 
charges fij. hj is the charge on grain j and <pj is the 
phase of the superfiuid order parameter. This model has 
a superfiuid ground state if the Josephson couplings dom- 
inate. With increasing charging energies it undergoes a 
quantum phase transition to an insulating ground state. 
In addition to Josephson junction arrays, the Hamilto- 
nian (TTJ describes a wide variety of other systems that 
undergo superfluid-insulator transitions. 

Within the strong-disorder approach [12J, the type of 
insulator depends on the symmetry properties of the off- 
set charge distribution. In contrast, these details were 
found unimportant at the critical point. In the follow- 
ing, we therefore focus on purely off-diagonal disorder, 
fij = 0. In this case, the Hamiltonian (fTJ) can be mapped 
onto a classical (1 + l)-dimensional XY model [3] 

H c \ = - ^ [Jj COs(<£j + l jT - 4> hT ) + Jj COSO^+I - jVr )] 
3,T 

(2) 

where j and r index the lattice sites in the space and 
time-like directions, respectively. The coupling constants 
Jj j T and Jj / T are determined by the parameters of the 
original Hamiltonian (fTJ with T being an effective "clas- 
sical" temperature, not equal to the real physical tem- 
perature which is zero. In the following, we fix J| and 
Jj and drive the XY model @ through the transition 
by tuning T. The interactions JJ and/or Jj are inde- 
pendent random variables drawn from probability distri- 
butions P S (J S ) and P*(J*). They depend on the space 
coordinate j only; this means the disorder is columnar 
(perfectly correlated in time direction). 

To determine the critical behavior of the classical XY 
model ©, we performed large-scale Monte-Carlo simula- 
tions using the efficient Wolff cluster algorithm 15j. We 
studied square lattices with linear sizes up to L = 3200 
and averaged the results over large numbers (200 to 3000, 
depending on L) of disorder realizations. Each sample 
was equilibrated using 200 to 400 Monte-Carlo sweeps, 
i.e., total spin flips per site. (The actual equilibration 
times both above and at the critical temperature T c did 
not exceed about 20 sweeps.) During the measurement 
period of 5000 to 30000 sweeps, we calculated observables 
such as specific heat, magnetization, susceptibility, spin- 
wave stiffness as well as correlation functions. In most 
simulations, we employed a uniform J? = 1 and drew the 
Jj from a binary probability distribution 



P*(J*) = cS( J 1 - r) + (1 - c)5{J l - 1) . 



(3) 



Here, c is the concentration of weak bonds which we 
fixed at c = 0.8. The disorder strength was tuned 
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FIG. 2. (Color online) Spatial correlation length £ s vs. tem- 
perature T for disorder strength r = 0.85 and system sizes 
L = 200 to 3200. The solid line is a fit to the KT form @. 
Inset: Luttinger parameter g at T c vs. system size L. 



by changing the value r of the weak bonds. In ad- 
dition to the clean case r = 1 (which corresponds to 
the pure superfhud-Mott insulator transition), we used 
r = 0.85,0.65,0.45,0.25, and 0.15. We also carried out 
test calculations with random J s . All simulations were 
performed on the Pegasus Cluster at Missouri S&T, using 
about 400,000 CPU hours 

We now turn to the results. To find T c for each dis- 
order strength r, we analyzed the behavior of the corre- 
lation length £ s (in the space-like direction indexed by 
j). It is calculated, as usual, from the second moment of 
the disorder-averaged correlation function. In the high- 
temperature phase but close to the transition, ^ s is ex- 
pected to follow the form 



6 = Aexp B{T~T C ) 



v-l/2 



(4) 



both in the clean KT universality class [lOj and in the 
strong-disorder scenario [l2j • A and B are non- universal 
constants. For all disorder strength, our data follow this 
prediction with high accuracy, see Fig. O for an example. 
We extract T c from fits of the data to ^ restricted to 
£, s > 10 to be in the critical region but £ s < L/10 to avoid 
finite-size effects. In the clean case (r = 1), we obtain 
T c = 0.8924 in excellent agreement with high-precision 
values in the literature [l6| [13] • 

In addition to the correlation length £ s in the space- 
like direction, we also studied the correlation length £ t in 
the time-like direction. We found £ t oc £ s for all disorder 
strengths which implies a dynamical exponent of z = 1. 

The order parameter susceptibility \ can b e analyzed 
analogously. In the high-temperature phase close to the 
transition, it is predicted to behave as 
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FIG. 3. (Color online) Susceptibility \ vs. temperature T for 
several disorder strengths. The maximum system sizes are at 
least L = 1500. The solid lines are fits to the KT form {SJ. 
The resulting estimates of T c are listed in the legend. 
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FIG. 4. (Color online) ln(x/ff) vs. ln(£ s ) for several disorder 
strengths and maximum system size L > 1500 (L — 500 for 
r = 0.15). The solid lines are linear fits; their slopes give —rj. 



Here, rj is the correlation function critical exponent and 
D = (2 — rj)B. Figure [3] shows that the data for all disor- 
der strengths r follow this prediction with high accuracy. 
The critical temperatures extracted from the correspond- 
ing fits are listed in the legend of the figure. Their values 
have small statistical errors ranging from about 3 x 10~ 4 
for the weak disorder cases to 2 x 10~ 3 for strong dis- 
order. The systematic errors due to corrections to the 
leading scaling form (J5J) are somewhat larger. We esti- 
mate them from the robustness of the fit against chang- 
ing the fit interval. This yields systematic errors ranging 
from about 5 x 1CP 3 for weak disorder to 2 x for 
strong disorder. Within these errors the critical temper- 
atures extracted from \ agree well with those from the 
correlation lengths. 

Equation ([SJ suggests a direct way to measure the ex- 
ponent 77: if one plots ln(x/£ 2 ) vs. ln(£ s ), the data should 
be on a straight line with slope —77. Figure [4] presents 
this analysis for different disorder strengths. In the clean 
case, r = 1, we find r/ = 0.243 in good agreement with the 
exact value 1/4 [10]. The weak-disorder curves (r = 0.85 
and 0.65) are parallel to the clean one within their sta- 
tistical errors. Fits in the range 20 < ^ s < L/10 give 
exponents 77 close to 1/4. In contrast, the strong-disorder 
curves (r = 0.45, 0.25, 0.15) are less steep, resulting in 
smaller 77. They are also noisier which leads to larger er- 
ror bars. All rj values are shown in Fig. [TJ They provide 
evidence for universal critical behavior (in the clean 2D 
XY universality class) for weak disorder but nonuniversal 
behavior for strong disorder. 

In addition to simulations in the high-temperature 
phase, we also studied the finite-size scaling properties 
of observables right at the critical temperature T c . Let 
us first consider the Luttinger parameter^ = n^/p s K. 



Under the quantum-to-classical mapping TJ], the com- 
pressibility k of the quantum rotor Hamiltonian ([1} maps 



onto the spin-wave stiffness pt in the time-like direction 
of the classical XY model ([2]). In our simulations, the 
Luttinger parameter is thus given by 



g = {-K/T)yfp^~ t 



(6) 



The stiffnesses p s and pt are not calculated by actually 
applying twisted boundary conditions during the sim- 
ulation but by using the relation given by Teitel and 
Jayaprakash [3] (for a derivation see, e.g., Ref. flit). 

Within KT theory, the Luttinger parameter close to 
the transition behaves as g(T) = g(T c ) + a(T c - T) 1 / 2 
where a is a constant and T < T c . Together with Q, 
this suggests the leading finite-size corrections to g at T c 
to take the form 



g(T c ,L) = g(T c ^) + b/\n(L) 



(7) 



where b is another constant. Calculating the Luttinger 
parameter at T c for different system sizes and extrapo- 
lating using ([7]) yields the infinite-system value g(T c , 00) 
[20I ] . We performed this analysis for all disorder strengths 
r and found that the g vs. l/ln(L) data indeed fall onto 
straight lines (the inset of Fig. [5] shows an example). The 
resulting extrapolated values are displayed in Fig. [TJ For 
weak disorder (r = 0.85 and 0.65), the Luttinger param- 
eters at T c agree with the clean value, g = 2, within 
their error bars (which are combinations of the statisti- 
cal Monte-Carlo error and the uncertainty in T c ). For 
stronger disorder (r = 0.45, 0.25, 0.15), g(T c , 00) takes 
larger, disorder-dependent values. 

Finally, we turn to the finite-size behavior of the sus- 
ceptibility at T c . According to finite-size scaling, the 
leading size-dependence should be of the form 



X (T C ,L)~L 2 -" 



(8) 
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behavior compatible with the clean 2D XY universality 
class, in agreement with Ref. 
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FIG. 5. (Color online) Susceptibility at T c plotted as \n{x/L 2 a ) 
vs. ln(L) for several disorder strengths. The solid lines are 
linear fits; their slopes give —r\ (values are shown in Fig. [TJ . 
The inset demonstrates the change of slope with increasing r. 



which provides another way to measure r/. Figure[5]shows 
plots of \n(x/L 2 ) vs. ln(L) for all disorder strengths r. 
For weak disorder (r = 0.85 and 0.65), the resulting val- 
ues of the exponent rj are again close to the clean value 
1/4. For larger disorder (r = 0.45, 0.25, and 0.15), we 
find disorder-dependent values that roughly agree with 
those extracted in the high-temperature phase (Fig. 0]). 

In summary, we used large-scale Monte-Carlo simu- 
lations to investigate the superfluid-insulator quantum 
phase transition of one-dimensional bosons with off- 
diagonal disorder. For weak disorder, our data provide 
evidence for a KT critical point in the universality class 
of the clean (1+1) dimensional classical XY model, with 
universal critical exponents rj = 1/4 and z = 1 as well as 
a universal value g = 2 of the critical Luttinger param 



eter. These results agree with the Harris criterion [21 1 
which predicts weak disorder to be an irrelevant pertur- 
bation at the clean KT transition. For stronger disor- 
der, the universality class of the transition changes. It 
is still of KT-type [£ s and \ follow (@J and ©] but the 
critical exponent r\ and the critical Luttinger parame- 
ter become disorder-dependent (non- universal) [22| . This 
agrees with the strong-disorder scenario [l2| ■ 

The important question of whether the boundary be- 
tween the weak and strong disorder regimes is sharp or 
just a crossover cannot be finally decided by means of our 
current numerical capabilities. The data in Fig. Q] would 
be compatible with both scenarios within their error bars. 

Earlier Monte- Carlo simulations [3] did not observe 
the strong-disorder regime. We believe that the binary 
disorder used in [l3} (equivalent to disorder in 3 s with 
c = 0.5 and r = 0.33 in our model) may not have been 
sufficiently strong. In particular, c = 0.5 is much less fa- 
vorable for the formation of rare regions than our c = 0.8. 
To test this hypothesis, we performed a few simulation 
using c = 0.5 and r — 0.33. They resulted in critical 



It is interesting to ask whether the different critical be- 
haviors in the weak and strong-disorder regimes are ac- 
companied by qualitative differences in the bulk phases. 
In particular, arc there two different insulating phases or 
are the weak and strong-disorder regimes continuously 
connected? A detailed analysis of the insulating phase(s) 
will also shed light on the mechanism that destroys the 
superfluid stiffness above T c . Is it due to the prolifer- 
ation of single quantum phase slips as at a clean KT 
transition or due to the formation of phase slip "dipoles" 
as suggested in Ref. [13]? Simulations to address these 
questions are under way. 

All our explicit results are for off-diagonal disorder and 
large commensurate filling. They do not directly apply to 
the generic dirty-boson problem with diagonal disorder 
considered in [j| [24 |. Note, however, that the critical 
behavior does not depend on the disorder type within 
the strong-disorder scenario Simulating the generic 
case would require a different approach (such as the link- 
current formulation 14() because the mapping onto a 
classical XY model is not valid for diagonal disorder. 

Finally, we turn to the experimental accessibility of 
the weak and strong-disorder regimes. Our results show 
that the transition between them occurs at a moderate 
disorder strengths. We therefore expect both regimes to 
be accessible in principle in experiments on systems such 
as ultracold atoms or Josephson junction chains (see also 
Ref. [H). 
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